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Abstract 

An  Eulerian,  sharp  interface,  Cartesian  grid  method  is  developed  for  the  numerical 
simulation  of  the  response  of  materials  to  impact,  shocks  and  detonations.  The  mass, 
momentum,  and  energy  equations  are  solved  along  with  evolution  equations  for 
deviatoric  stresses  and  equivalent  plastic  strain.  These  equations  are  cast  in  Eulerian 
conservation  law  form.  The  Mie-Gruneisen  equation  of  state  is  used  to  obtain  pressure 
and  the  material  is  modeled  as  a  Johnson-Cook  solid.  The  ENO  scheme  is  employed  to 
capture  shocks  in  combination  with  a  hybrid  particle  level  set  technique  to  evolve  sharp 
immersed  boundaries.  The  numerical  technique  is  able  to  handle  collisions  between 
multiple  materials  and  can  accurately  compute  the  dynamics  of  the  immersed 
boundaries.  Results  of  calculations  for  axisymmetric  Taylor  bar  impact  and  penetration 
of  a  Tungsten  rod  into  steel  plate  show  good  agreement  with  moving  finite  element 
solutions  and  experimental  results.  Qualitative  agreement  with  theory  is  shown  for  the 
void  collapse  phenomenon  in  an  impacted  material  containing  a  spherical  void. 

1.  Introduction 

An  Eulerian  methodology  is  presented  for  computing  a  range  of  problems  that  can  be 
classified  as  high-speed  multimaterial  interactions(Zukas,  1982,  Meyers,  1994).  Such 
interactions  can  arise  in  phenomena  such  as  munition-target  interactions,  geological 
impact  dynamics,  shock-processing  of  powders,  formation  of  shaped  charges  upon 
detonation  and  their  subsequent  interaction  with  targets,  and  the  impact-induced 
detonation  of  porous  high-explosives. 

The  fundamental  challenges  to  a  simulation  capability  designed  to  solve  problems  in  the 
physical  phenomena  listed  above  arise  from  the  presence  of  nonlinear  wave  propagation 
and  the  large  deformations  suffered  by  the  interacting  materials.  Therefore,  to  simulate 
the  physics  of  multimaterial  interactions,  the  numerical  approach  should  be  able  to  track 
the  propagating  boundaries  and  shock  waves  simultaneously  and  accurately. 
Traditionally,  methods  that  have  been  used  to  solve  such  problems  have  been  termed 
hydrocodes.  The  broad  range  of  available  hydrocodes  has  been  reviewed  by  Anderson 
(1987)  and  Benson  (1992).  The  wide  variety  of  available  methods  indicates  both  that 
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computational  research  in  this  area  has  been  a  robust  enterprise  and  that  the  existing 
methods  all  carry  limitations  that  have  spurred  continuing  development  efforts. 

Hydrocodes  treat  the  moving  material  boundaries  by  either  allowing  the  boundaries  to 
flow  through  a  fixed  mesh  while  computing  the  flow  field  on  the  fixed  mesh,  or  by 
allowing  the  mesh  to  follow  the  material  points  in  the  deforming  materials.  An 
intermediate  approach,  ALE  (Arbitrary  Lagrangian  Eulerian),  allows  the  mesh  to  move 
so  as  to  conform  to  the  contours  of  the  deforming  object,  but  the  mesh  is  not  necessarily 
attached  to  the  material  points.  Lagrangian  and  Arbitrary  Lagrangian  Eulerian  (ALE) 
methods  (Liu  et  al.,  1986)  for  the  simulation  of  problems  with  severe  material 
deformation  have  been  applied  extensively  in  the  solid  mechanics  community.  For 
example,  Camacho  and  Ortiz  (1996,  1997)  have  performed  Lagrangian  finite  element 
calculations  for  impact  and  deformation  of  brittle  materials  (Camacho  and  Ortiz  1996), 
and  ductile  penetration  (Camacho  and  Ortiz  1997).  Their  approach  is  based  on  adaptive 
meshing,  explicit  contact/friction  algorithm,  and  rate  dependent  plasticity.  In  moving 
mesh  methods,  considerable  complexity  is  enjoined  by  the  need  for  mesh  management, 
i.e.  in  maintaining  an  adequately  refined  mesh  with  good  mesh  quality.  For  very  severe 
deformations,  meshless  methods  (Duarte  and  Oden,  1996,  Johnson  et  al.,  1996,  Liu  et  al., 
2000,  Belytschko  et  al.,  2000),  or  a  combination  of  finite  element  methods  with 
embedded  boundary  tracking  and  local  enrichment  (Dolbow  et  al.,  2000,  Moes  et  al., 
1999,  Sukumar  et  al.,  2000),  have  emerged  as  attractive  alternatives  in  recent  years.  In 
these  methods,  one  either  entirely  dispenses  with  a  mesh  or  the  mesh  does  not  distort  as 
the  embedded  boundary  (such  as  a  crack)  propagates  through  the  mesh. 

Eulerian  methods  have  been  applied  to  study  material  deformation  by  some  researchers 
by  adapting  techniques  from  computational  fluid  dynamics.  For  example,  Trangenstein 
(1990,  1994,  1995),  and  Trangenstein  and  Pember  (1991)  have  adopted  Godunov’s 
method  to  handle  multimaterial  impact  as  a  Riemann  type  problem  with  second-order 
accuracy.  Benson  and  coworkers  (Benson,  1995,  Cooper  et  al.,  2000)  have  applied 
Eulerian  methods  to  study  the  collapse  of  voids  and  shock-induced  compaction  in 
materials  under  impact  loading.  The  methods  presented  by  Benson  and  coworkers, 
although  based  on  an  Eulerian  fixed  mesh  setting  are  of  the  Lagrangian-plus-remap  type, 
where  the  material  deformation  calculations  are  split  into  two  steps.  First  the  material  is 
evolved  by  a  Lagrangian  step  which  deforms  nodes  to  new  positions  and  then  the  field  is 
mapped  back  to  the  fixed  Eulerian  mesh  and  the  new  interfaces  reconstructed  using 
Young’s  reconstruction.  This  approach  has  been  used  to  good  effect  in  the  solution  of 
mesoscale  response  of  energetic  materials  in  shock  compression  (Menikoff  and  Kober, 
1999)  and  void  collapse  (Cooper  et  al.,  2000). 

Shock-capturing  methods  that  were  developed  for  gas  dynamics  have  been  extended  to 
condensed  media  for  application  to  high  velocity  (in  liquids)  or  high  strain  rate  (in  solids) 
problems  where  nonlinear  wave  propagation  phenomena  are  important.  Glaister  (1988) 
and  Arienti  et  al.  (1999)  employ  the  Roe  scheme  in  an  approximate  Riemann  solver  to 
capture  shocks.  While  the  former  work  is  restricted  to  gases  and  1 -dimension  with  a 
general  convex  equation  of  state,  the  latter  deals  with  solid  materials  with  the  Mie- 
Griineisen  equation  of  state  for  the  pressure,  but  they  solve  the  Euler  equations  for  the 
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flow  of  the  condensed  material,  i.e.  the  strength  of  the  material  is  not  considered.  Arienti 
et  al.  (1999)  have  also  investigated  2-dimensional  problems  in  that  setting.  Following 
Dukowicz  (1985),  Miller  and  Puckett  (1996)  presented  an  approximate  Riemann  solver 
for  multimaterials  for  general  equations  of  state  where  material  interfaces  can  lie  within 
cells.  They  treated  the  multiple  materials  as  a  mixture  within  each  cell  (i.e.  volume 
fractions)  without  resorting  to  the  Lagrangian-plus-remap  approach.  Material  strength 
was  not  considered.  The  discrete  Riemann  solver  for  their  formulation  was  fairly 
challenging  to  develop,  particularly  at  the  faces  of  the  mixed  cells.  A  simpler  approach 
is  the  Ghost  Fluid  Method  due  to  Fedkiw  and  coworkers  (1999a).  In  this  method  the 
interface  is  treated  as  a  sharp  entity  that  resides  on  the  fixed  mesh  and  appropriate 
boundary  conditions  at  the  interface  are  applied  by  extrapolating  the  field  to  an  extended 
“ghost”  material.  This  approach  leads  to  a  local  reduction  in  order  of  accuracy  at  the 
computational  points  adjoining  the  immersed  interfaces.  However,  since  such  points  are 
few  in  number,  the  overall  accuracy  is  still  maintained  at  the  high-order.  Fedkiw  et  al. 
(1999a)  have  applied  the  ENO  schemes  to  study  the  propagation  of  shocks  in  media 
where  the  pressure  is  governed  by  a  variety  of  equation  of  states  (gases  and  liquids). 
Their  results  show  that  the  ENO  scheme  can  accurately  handle  the  shock  formation  in 
such  systems. 

In  a  recently  published  paper  (Udaykumar  et  al.,  2003)  we  presented  an  Eulerian 
methodology  to  simulate  high-speed  impact  of  materials.  The  ENO  scheme  was  used 
along  with  explicit  interface  tracking  and  a  Prandtl-Reuss  material  model  to  describe 
elasto-plastic  deformation.  Validation  exercises  were  carried  out  for  one-dimensional 
hydrodynamic  and  elasto-plastic  impact  situations.  These  ID  results  showed  good 
agreement  with  exact  solutions  and  with  a  Lagrangian  technique  using  a  moving  mesh. 
Although  a  third-order  ENO  scheme  was  used  the  fact  that  the  material  was  only  weakly 
compressible  rendered  the  captured  discontinuities  to  be  somewhat  more  smeared  than 
for  strong  hydrodynamic  shocks.  Sharpening  of  the  shock  by  use  of  different  limiters 
along  with  the  ENO  scheme  was  studied.  Based  on  this  study,  the  Local-Lax  Friedrichs 
ENO  scheme  with  the  minmod  limiter  was  found  to  be  the  most  robust  and  captured  the 
physical  features  of  elasto-plastic  waves  with  optimal  resolution.  The  present  paper 
advances  the  methodology  by  employing  the  hybrid  particle  level  set  method  to  track 
boundaries  and  by  incorporating  rate-dependent  plasticity  (through  the  Johnson-Cook 
model)  to  better  represent  the  dynamics.  In  addition,  substantial  improvements  and 
simplifications  in  application  of  interface  conditions  are  derived  from  the  level-set 
representation. 

The  characteristics  of  the  present  numerical  method  that  make  it  attractive  relative  to 
other  methods  employed  in  hydrocodes  (including  our  previous  work  presented  in 
Udaykumar  et  al.  (  2003))  for  high-speed  multimaterial  flows  are: 

1.  The  equations  governing  the  material  deformation  are  solved  in  an  Eulerian 
setting  on  a  fixed  Cartesian  mesh.  Well-developed  high-accuracy  shock  capturing 
schemes  are  easily  applied  to  compute  the  nonlinear  wave-propagation 
phenomena  in  this  framework.  Here  the  ENO  scheme  is  used  for  all  calculations. 
Addition  of  problem-dependent  shock  viscosity  is  not  called  for  since  adequate 
dissipation  at  discontinuities  is  built  into  the  scheme. 
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2.  The  mesh  remains  fixed  while  the  material  flows  through  the  mesh.  Thus,  issues 
such  as  mesh  deformation,  entangling,  catastrophic  mesh  distortion  in  regions 
that  have  changed  phase  from  solid  to  liquid  and  thus  have  lost  strength,  do  not 
arise  within  the  current  fixed-grid  approach.  The  materials  can  fragment  or 
collide  and/or  merge  without  affecting  the  flow  calculations. 

3.  The  interfaces  are  tracked  in  a  sharp  fashion.  They  are  not  smeared  over  the  mesh 
as  in  traditional  Eulerian  methods.  Thus  materials  can  approach  each  other 
without  mixing  and  a  mixture  formulation  is  not  required  in  treating  cells  with 
multiple  interfaces  or  in  cells  that  are  only  partially  filled.  The  exact  interface 
location  is  known  at  all  times  due  to  the  level  set  representation.  Boundary 
conditions  and  jump  conditions  can  be  applied  at  the  sharp  interfaces  at  both  free 
surfaces  and  material-material  interfaces. 

4.  A  particle-level  set  method  (Enright  et  al.,  2001)  is  used.  This  technique  is  shown 
to  maintain  sharp  comers  of  objects  without  excessive  smoothing  due  to  the 
inherent  entropy  fix  in  the  advection  scheme  used  to  evolve  the  narrow-band  grid- 
based  level  set.  Thus,  spurious  mass  loss  effects  due  to  stretching  and  shearing  of 
interfaces  that  plague  all  Eulerian  interface  tracking  schemes  are  avoided  in  the 
present  method.  No  difficulty  arises  in  treating  multiple  boundaries.  These  are 
simply  evolved  as  different  level  set  functions.  The  interfaces  can  undergo 
topological  changes  without  occasioning  difficulties  for  the  flow  solver. 

5.  Extension  to  3-dimensions  should  be  straightforward.  The  numerical  schemes  for 
the  governing  equations  are  obtained  by  field-by-field  decomposition  along  each 
dimension  and  therefore  addition  of  the  third  dimension  will  only  necessitate 
discretization  of  the  equations  in  that  direction.  Furthermore,  the  level  set 
formulation  is  also  easily  extended  to  3D,  thus  allowing  for  tracking  of  three- 
dimensional  objects  as  a  straightforward  extension  of  the  technique  presented  in 
this  work. 

The  method  is  benchmarked  in  this  work  by  comparing  with  solutions  using  moving 
FEM  techniques  and  experimental  data. 

2.  Formulation  of  the  Problem 

2.1  Constitutive  Relations 

The  equations  governing  the  material  deformation  appropriate  for  high  strain  rate 
applications  can  be  formulated  by  assuming  that  the  volumetric  or  dilatational  response 
is  governed  by  an  equation  of  state  while  the  deviatoric  response  obeys  a  conventional 
flow  theory  of  plasticity.  The  system  of  equations  describing  the  material  deformation  in 
Lagrangian  form  can  be  written  as  follows.  The  stress  in  the  material  is  expressed  as  the 
sum  of  the  dilatational  and  deviatoric  parts: 

<r..  =sy  -pSy  (1) 

Here  aiy  is  the  Cauchy  stress  tensor,  stj  its  deviatoric  part,  and  p  the  hydrostatic  pressure 

taken  to  be  positive  in  compression. 

The  rate  of  change  of  deviatoric  stress  is  given  by: 

1,=2G{D9-D')  (2) 
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Sy  =  s(j  +  SfcQy  —  ^iicskj  (3) 

v 

where  G  is  the  shear  modulus,  sy  the  Jaumann  derivative  and  sy  the  total  derivative  of 
Sy,  Dy  the  deviatoric  strain-rate,  Dp  the  plastic  strain-rate,  and  Qy  the  spin  tensor.  The 

Jaumann  derivative  is  used  to  ensure  material  frame  indifference  with  respect  to  rotation. 
The  spin  tensor  and  the  strain-rate  tensor  are  given  by: 

<4) 

Dy-\iy,j^hi)  (5) 

where  Dy  =  D?  +  Djj  is  the  strain-rate  tensor  composed  of  the  elastic  strain  rate  Dy  and 

Qy 

the  plastic  strain  rate  D?  and  v*  is  the  ith  velocity  component  and  v,  y  =  . 

The  deviatoric  strain-rate  is: 

D,=D,-±DaS,  (6) 

The  plastic  strain-rate  follows  the  relationships: 

tr(D?)  =  0  (7) 

O'  =  A N,  (8) 

s ..  .  . 

where  Ntj  =  l  -  is  the  unit  outward  normal  to  the  yield  surface  and  A  a  positive 

slSUSkl 

parameter  called  the  consistency  parameter. 

The  effective  stress  (Se)  and  effective  plastic  strain  (ep )  are  given  by: 


- 1 <Kv>) 

(9) 

{Pf  =|  hd;d;) 

(10) 

The  evolution  of  the  temperature  due  to  heat  conduction  and  thermal  energy  produced  by 
work  done  during  elasto-plastic  deformation,  is  written  as: 

pCT  =  kV2T-a{M  +  2ii)T0eekk+pWp  (11) 

where  T  is  the  temperature,  C  the  specific  heat,  k  the  thermal  conductivity,  a  the 
thermal  expansion  coefficient,  A.  and  /u  are  the  Lame  constants,  p  the  Taylor-Quinney 
parameter  implies  the  fraction  of  mechanical  power  convert  to  thermal  power  (Taylor  and 
Quinney,  1934),  and  is  taken  as  0.9,  and  the  stress  power  due  to  plastic  work, 
Wp  =spSe.  For  the  applications  considered  in  this  work,  the  conduction  and  elastic 
work  terms  are  small  compared  to  the  plastic  work  term. 

The  material  models  used  in  this  work  are  rate-independent  (Prandtl-Reuss  material)  and 
shear  rate  and  temperature  dependent  (Johnson-Cook,  1985)  models,  given  (respectively) 
in  general  form  as: 

<jy=A  +  B(£PY  (12) 
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T-T 

Here  A,B,C,n,m,s0  are  model  constants,  and  8  = - — ,  where  T0  and  Tm  are 

reference  room  temperature  and  melting  temperature,  respectively.  In  the  Johnson-Cook 
model,  equation  (13),  the  flow  stress,  cr  v ,  increases  with  an  increase  in  the  effective 

plastic  strain  and  the  effective  plastic  strain-rate,  and  decreases  with  an  increase  in 
temperature.  The  yield  stress  in  fact,  goes  to  zero  as  the  temperature  approaches  the 
melting  temperature  of  the  material. 


2.2  Governing  Equations 


In  the  Eulerian  setting,  since  material  points  are  not  followed,  the  above  constitutive 
equations  need  to  be  combined  with  calculations  of  the  flow  of  material.  The  equations 
governing  axisymmetric  deformation  and  flow  of  the  material  can  be  written  as  a  system 
of  conservation  laws  using  the  primitive  variables  as: 

(.4) 

dt  dx  dy 

where  the  vector  of  conserved  variables  ( Q ),  and  the  convective  flux  vectors  in  the 


The  above  equations,  written  in  conservative  form  include  the  evolution  of  the  deviatoric 
stresses  according  to  Eq.  (2).  The  dilatational  part  of  the  stress  (pressure)  is  obtained 
from  an  equation  of  state.  The  source  vector  that  appears  in  Eq.  (14)  is  of  the  form: 
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(16b) 

(16c) 
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where  E  =  — 
3 


lfdu  u  dv  ')  ^^-Qorl  for  the  2D  or  2D-axisymmetric  formulations 


KS  W  V  F 

- ! - 1 - 

dx  x  dy  j 


respectively. 

The  equation  for  evolution  of  the  plastic  strain  and  temperature  fields  are  also  solved 
deduce  the  thermo-mechanical  effects  of  the  multimaterial  interactions.  These  are  wri 
as: 

ds” 


dt-  +  V-Ve»=^-tr(D'Dp 
-  +  F-VT  =  \{i W2T  - a(3A  +  2 +  PWP ) 

pc 


to 

are  written 
(17) 


K  +  v-VT  =  —  (yfcV27’-a(3A  +  2//)r04  +  pWp)  (18) 
dt  pC 

The  eigenvalues  of  the  equation  system  (14)  were  found  to  be  real  (Udaykumar  et  al., 
2003)  for  the  range  of  parameters,  i.e.  material  properties  and  impact  velocities,  of 
interest  in  this  work  For  ID  formulation,  the  eigenvalues  were  shown  to  be  (Vanden, 
1998): 

A,  =A2  =A3  =u  (19a) 
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As  =c  —  u  (19c) 

The  sound  speed  depends  on  the  particular  equation  of  state  chosen  for  the  pressure  and 
is  given  below  for  the  case  of  the  Mie-Griineisen  equation  of  state. 

The  above  system  of  equations  are  first  solved  using  the  Essentially  Non-Oscillatory 
(ENO)  numerical  scheme  (Shu  and  Osher,  1988,  1989)  as  will  be  described  in  a 
subsequent  section.  The  radial  return  algorithm  is  then  applied  to  bring  the  stress  state 
back  to  the  yield  surface  if  the  predicted  von  Mises  stress  falls  outside  the  admissible 
states.  Equations  (17)  and  (18)  for  the  effective  plastic  strain  and  temperature  field  are 
solved  separately  using  a  simple  second-order  upwind  scheme. 

2.3  Mie-Griineisen  e.o.s.  for  Solids 


It  is  necessary  to  relate  pressure,  specific  volume,  and  internal  energy  through  an 
equation  of  state  in  order  to  close  the  system  of  equations  above.  The  overall  pressure  in 
a  solid  may  be  expressed  as  a  sum  of  two  terms,  the  first  due  to  the  thermal  excitation, 
and  the  second  due  to  the  pressure  resulting  from  the  attractive/repulsive  forces  in  the 
lattice  of  atoms.  For  the  present  work,  due  to  the  high  strain-rate,  large  deformation 
problem  of  concern,  we  utilize  the  Mie-Griineisen  e.o.s.  If  ec  and  pc  denote  respectively 
the  “cold”  (at OAT)  energy  and  pressure,  the  incomplete,  temperature-independent 
formulation  of  the  Mie-Griineisen  e.o.s.  is: 


(e-e(V)) 

P(e,  V)  ~  r(K) - ^ —  +  PC(V)  =  r-  +  f(V) 

where,  by  definition, 


P  2 
l 

V  =  — 
p 

and  the  Griineisen  parameter  is  defined  as: 


(20) 

(21) 

(22) 


r  =  v\ 


Vo 


(23) 


'V  p 

where  po  is  the  density  of  the  unstressed  material,  Co  and  s  are  coefficients  that  relate  the 
shock  speed  Us  and  the  particle  velocity  up .  Experiments  on  solids  provide  a  relation 

between  Us  and  up  (Olinger  et  al.,  1975).  A  first  approximation  consists  of  a  linear 

relation  (generally  applicable  for  strong  shocks): 

U5=c0+sup  (24) 

where  c0  is  the  bulk  sound  speed  for  the  material  at  rest,  and  5  is  related  to  the 
isentropic  pressure  derivative  of  the  isentropic  bulk  modulus.  In  general,  c0  and  s  are 
obtained  from  experiments. 


The  expression  for  the  sound  speed  is  given  by: 


c2  =f^l  +JL($pX  =Te  +  f\V)+r?-  (25) 

\dp)e  p  \de)p  p 

For  low  energy  and  low  pressure  (as  in  a  rarefaction),  round-off  and  approximation  errors 
may  result  in  negative  values  of  c2,  since  /'  changes  sign.  To  correct  this,  and 
following  the  approach  by  Miller  and  Puckett  (1996),  Arienti  et.  al.  (1999)  prolong  the 
e.o.s.  with  a  pseudo  elastic-solid  e.o.s.  The  final  expression  for  f(V) ,  to  accommodate 


for  negative  pressure  (tension)  and  preserve  positive  sound  speed-squared,  is  written  as: 

t r%[1-£(F»_K)l  if  v~v° 

.  2V  J 


fiV)  = 


if  V<V0 


I__L 

V  vn 


if  V>V0 


where  q>~\- 


2.4  Radial  Return  Algorithm 

When  a  material  deforms  plastically,  the  set  of  stresses  marking  the  transition  boundary 
between  the  admissible  and  the  inadmissible  stress  states  is  defined  by  the  yield  surface. 
The  stress  must  be  constrained  to  always  fall  either  within  or  on  the  yield  surface. 
Typically,  numerical  solutions  of  the  plasticity  equations  tentatively  assume  that  the 
material  response  for  the  entire  time  increment  is  purely  elastic.  If  the  resulting 
“predicted”  or  “trial”  updated  stress  is  found  to  fall  outside  the  yield  surface,  then  the 
numerical  algorithm  recognizes  that  the  tentative  assumption  of  elasticity  must  have  been 
wrong.  Classical  return  algorithms  assert  that  the  correct  updated  stress  can  be  obtained 
by  projecting  the  inadmissible  trial  stress  back  to  the  yield  surface.  There  are  various 
methods  by  which  such  return  of  the  stress  to  the  yield  surface  can  be  effected  (Montan 
and  Boija,  2002,  Bilotta,  2001,  van  Le  and  Saxce,  2000).  Here,  we  adopt  an  algorithm 
due  to  Ponthot  (1998). 

We  assume  the  existence  of  a  general  yield  function  /  given,  for  a  J2  von  Mises  material 
with  isotropic  hardening,  by: 

f(s,,S.)  =  S.-o,  (27) 

where  crv  is  the  current  yield  stress. 

The  admissible  stress  states  are  constrained  to  remain  on  or  within  the  domain  defined  by 
this  yield  function,  i.e.  we  require  that  /  <  0 .  The  hardening  law  is  given  by: 
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where  h  is  called  the  hardening  coefficient  and  corresponds  to  the  slope  of  the  effective 
stress  versus  effective  plastic  strain  curve  under  uniaxial  loading  conditions.  Substitution 


^  [2 

of  Eq.  (8)  into  Eq.  (10)  results  in  s  p  =  J-A .  Eq.  (28)  can  also  be  written  as: 

1  3 


dv  =  hP  (29) 

which  clearly  indicates  the  role  of  the  hardening  coefficient.  Time  integration  of  the 
deviatoric  stress  equation  for  elasto-plastic  deformations,  i.e.: 

*,  +•*»£>*  =  20(5 -D')  (30) 

is  split  into  two  parts,  the  elastic  predictor  which  yields  a  trial  stress  assuming  purely 
elastic  deformation  of  the  material: 


Sij.tr  +  Sik,tr^kj  &ikS kj,tr  ~  ^GD  (31a) 

and  the  plastic  corrector  to  bring  the  computed  trial  stress  back  to  the  yield  surface: 

Sy.cor  =  -2GD'  =  -2GAN,  (31b) 

where  the  subscript  tr  is  for  trial,  and  cor  is  for  correction.  Note  that  the  update  of  the 
deviatoric  stress  in  Eq.  (31a)  is  reflected  in  the  source  terms  in  Eq.  (16b-d).  Following  the 
update  using  the  system  of  equations  in  Eq.  (14)  the  correction  of  the  deviatoric  stresses 
through  the  radial  return  is  performed  using  Eq.  (31b).  The  unit  normal  to  the  yield 
surface  is  approximated  as: 


The  return  is  effected  in  a  direction  normal  to  the  yield  surface  as  follows: 

Sy.cor  ~  Sy.tr  ~  2G4N jjjr 


(32) 

(33) 


+J7  (34) 

'l 

where  the  unknown  scalar  parameter  £  =  Ja dt  where  to  and  ti  are  the  beginning  and  end 

of  the  time  interval  of  integration.  This  parameter  is  determined  by  the  enforcement  of 
the  consistency  condition,  /  =  0 ,  at  time  t  =  tx ,  i.e.  we  require  that: 

/-■Jf  k,  -2G^s,rIv -IGCN^-o'.  =0  (35) 

Eq.  (35)  can  be  solved  for  £  by  first  discretizing  Eq.  (28)  in  time  as: 


where  superscripts  0  and  1  denote  the  values  at  to  and  ti  respectively.  Substituting  for  cr‘ 
using  Eq.  (36)  into  Eq.  (35)  and  solving  for£ ,  one  obtains  the  final  expression: 
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Thus,  once  £is  obtained  the  correction  of  the  predicted  deviatoric  stresses  is  performed 
using  Eq.  (31b)  and  the  consistency  condition  is  enforced. 

3.  Numerical  Method 

3.1.  Local  Lax  Friedrichs  Essentially  Non-Oscillatory  (LLF-ENO)  Scheme 

To  solve  the  hyperbolic  system  of  equations  in  1-  and  2-dimensions,  the  ENO  shock¬ 
capturing  scheme  (Shu  and  Osher,  1988,  1989)  is  used.  The  Convex  ENO  scheme  due  to 
Liu  and  Osher  (1998)  is  implemented,  to  enable  the  oscillation-free  solution  of  the  2- 
dimensional  equations  without  field-by-field  decomposition  in  the  presence  of  large 
gradients.  The  discretization  has  been  described  in  detail  in  a  previous  paper  (Udaykumar 
et  al.,  2003)  and  is  only  briefly  outlined  below. 

Consider  the  governing  equation  for  one-dimensional  transport: 


^+af(e)=s(e) 

dt  8x 

(38) 

Let 

f-KQ) 

(39) 

where 

Q 

+ 

i 

1 

II 

(40) 

*'-*w 

Fe  and  Fw  are  the  fluxes  at  the  east  and  west  faces  shown,  and  xe  and  xw  are  the  locations 
of  the  east  and  west  faces  respectively,  as  shown  in  Figure  1(a).  D  is  an  appropriate 
discrete  operator  for  the  source  terms.  In  the  current  work,  the  source  terms  are 
discretized  using  a  2nd-order  central  difference  scheme.  This  was  found  to  be  robust  for 
the  calculations  performed.  However,  it  may  be  necessary  in  future  work  to  develop  a 
more  sophisticated  differencing  procedure  for  the  source  terms  as  well. 

The  three-step  third-order  in  time  Runge-Kutta  scheme  is  used  in  this  work  and  takes  the 
form  (Shu  and  Osher,  1988,1989): 

Qm  =  Q(n)  +  AtL(Qin)) 

Qm  =  +!Q(n))  +  ^tL(Q({))  (41) 

Q(n+\)  _  !(2g(2)  +Q(n))  +  -AtL(Q(2)) 

3  3 

The  spatial  order  of  accuracy  of  the  ENO  formulation  used  to  solve  Eq.  (38)  is 
determined  by  the  interpolation  practices  used  to  evaluate  the  fluxes  at  the  faces  e  and  w, 
i.e.  in  obtaining  Fe  and  Fw  in  Eq.  (40).  Due  to  the  presence  of  immersed  boundaries,  as 
illustrated  in  Figure  1(b),  care  must  be  taken  in  computing  these  fluxes.  The  flux 
evaluations  for  the  ENO  formulation  comes  from  derivatives  of  an  interpolating  function 
H(x)  as  follows: 

<42> 


11 


The  derivatives  are  evaluated  from  divided  differences  and  the  flux  evaluation  is 
performed  as  follows: 

Fe  =  f;  +  F;  =  H'+  (xe)  +  H-  (xj  (43) 

The  superscripts  (+)  and  (-)  indicate  the  positive  and  negative  direction  fluxes  at  the  face 
e  under  consideration  as  illustrated  in  Figure  1(a).  The  derivatives  H'  are  obtained  by 
polynomial  reconstruction  of  the  fields  as  in  the  standard  ENO  implementations.  The 
present  formulation  based  on  the  Convex  ENO  scheme  proposed  by  Liu  and  Osher 
(1998)  chooses  the  divided  difference  value  “closest”  to  the  previous  order  flux  chosen. 
The  scheme  reduces  to  low-order  automatically  at  discontinuities,  while  maintaining 
higher-order  in  smooth  regions.  The  first  divided  difference  is  obtained  as  follows: 

H'+(xe)  =  H+  x  ,,xe  =^(f(q[Xj])  +  a  ,?[x;])  (44a) 

L  J~i  J  2  J+2 

and 

H"(xe)  =  H-  xe,x  3  =Uf(q[xj+,])-a  tq[xj+ ,])  (44b) 

j+2 J  2  '+2 

where  ccj+V2  is  the  characteristic  speed  evaluated  at  the  cell  face  location  xJ+u2 .  This  is 
evaluated  as  the  maximum  eigenvalue  of  the  set  in  Eq.  (38)  at  the  cell  face. 

For  points  that  are  adjacent  to  the  immersed  interface  such  as  j  in  Figure  1(b),  the  flux 
evaluations  need  to  be  modified.  Here,  the  east  face  is  still  taken  at  the  grid  cell  face  and 
for  point  j: 

H'+(xe)  =  H[x  ,,xj  =  \(f(q[xj])  +  a  ^[xj])  (45a) 

J-2  2  J+~2 

H  (Xg )  —  H\Xe,X  3  ]  —  “ (/(<7 phantom  )  —  ®  ,+± Q phantom  )  (45b) 

where  qghost  is  the  ghost  value  of  the  convected  scalar  variable  q  (see  Fig.  1(c)).  This 
value  needs  to  be  obtained  while  satisfying  appropriate  boundary  conditions  on  the 
immersed  interface  as  described  below.  This  type  of  interfacial  flux  treatment  of  course 
reduces  the  order  of  accuracy  at  the  immersed  boundaries  by  one  order.  However,  the 
high-order  scheme  is  retained  in  the  bulk  of  the  computational  domain.  Similar 
considerations  apply  in  the  Ghost  Fluid  Method  (Fedkiw  et  al.,  1999a)  for  multifluid 
interactions,  and  flux  separation  method  (Tran  et  al.,  1999)  for  the  treatment  of  shock  and 
contact  discontinuities. 

3.2  Boundary  Conditions 

To  evaluate  the  fluxes  in  the  discrete  form  as  in  Eq.  (40),  appropriate  boundary 
conditions  need  to  be  applied  at  the  interface  location.  As  can  be  seen  in  Figure  2(a), 
when  two  interfaces  collide  a  portion  of  an  interface  in  the  material-material  contact 
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region  can  have  an  immediately  adjoining  portion  that  falls  in  the  material-void  region  of 
the  interface.  At  each  instant  of  the  interface  deformation,  the  ghost  points  (defined  as 
points  which  are  immediately  outside  a  given  body  as  indicated  in  Figure  2(b)) 
associated  with  each  portion  of  the  interface  have  to  be  classified  as  collision  (material- 
material  or  M-M)  or  free  (material-void  or  M-V)  ghost  points  and  the  appropriate  b.c.’s 
applied  at  the  interfaces  by  supplying  values  to  the  ghost  points.  The  procedure  for 
determining  how  a  ghost  point  of  an  interface  is  to  be  classified  as  collision  or  free  ghost 
point  is  discussed  in  section  4.3. 

At  the  outset,  all  the  dependent  variables  are  extrapolated  from  the  interior  of  a  material 
to  the  ghost  points.  Thereafter,  depending  on  the  boundary  conditions,  these  extrapolated 
values  are  replaced  by  values  that  satisfy  the  boundary  conditions  at  the  boundary  of  the 
material  in  question.  Therefore,  before  discussing  the  boundary  conditions,  it  is  useful  to 
outline  the  procedure  for  extrapolating  values  from  the  interior  of  a  given  material  to 
ghost  points  lying  across  the  boundary  of  the  material.  In  order  to  extrapolate  values 
from  an  object  to  the  ghost  point  P,  as  illustrated  in  Figure  3,  the  ghost  point  at  P  is  first 
reflected  (along  the  local  normal  obtained  from  the  level  set  field)  across  the  material 
interface.  The  values  of  all  the  dependent  variables  at  the  mirror  point  Ij  are  then 
determined  by  interpolation  from  the  surrounding  nodes.  Care  is  exercised  to  ensure  that 
there  are  four  immediate  surrounding  nodes  lie  in  the  same  phase  in  which  point  Ii 
resides.  If  this  is  the  case,  then  bilinear  interpolation  procedure  can  be  carried  out  in  a 
straightforward  way  to  obtain  values  for  point  Ii .  A  second  point,  Point  I2  in  Figure  3  is 

also  placed  at  a  distance  dl  =  sJax2  +  Ay2  from  point  Ii  (Figure  3)  along  the  normal  from 
ghost  point  P  to  the  interface.  Again  bilinear  interpolation  is  employed  to  obtain  values 
for  point  I2  if  there  exist  four  immediate  neighboring  nodes  in  the  same  phase.  If  all  four 
of  the  surrounding  points  are  not  in  the  same  material  for  either  point  Ii  or  I2,  then 
distance-weighted  interpolation  from  the  surrounding  grid  points  that  lie  in  the  same 
phase  is  performed  to  obtain  values  for  that  point.  Using  values  at  points  Ii  and  I2,  linear 
extrapolation  is  then  performed  to  obtain  values  for  the  ghost  point  P.  Initially  all  the 
dependent  variables  are  extrapolated  to  the  ghost  points.  These  values  are  then  replaced 
as  described  below  depending  on  the  type  of  boundary  condition  to  be  applied  on  the 
interface  adjoining  the  ghost  point. 

Type  1:  Material-material  interface 

Boundary  conditions  on  portions  of  the  interface  which  qualify  as  material-material 
interfaces  are  developed  based  on  the  physically  imposed  interface  conditions.  At  a 
material-material  (or  collision)  interface  we  enforce  continuous  material  point  velocities 
normal  to  the  interface  for  the  two  materials  and  the  continuity  of  normal  traction  and 
temperature,  whereas  the  tangential  traction  component  may  remain  discontinuous.  The 
procedure  for  applying  the  interface  conditions  is  illustrated  in  Figure  4(a),  where  the 
subscript  A  denotes  object  A  and  subscript  B  denotes  object  B.  We  want  to  determine 
the  variables  at  point  PA  (ghost  point  of  A,  which  lies  inside  object  B)  by  replacing  the 
extrapolated  values  with  ones  that  apply  the  required  boundary  conditions  at  the  interface 
location  O  shown  in  Figure  4(a).  First,  the  unit  normal  vector  is  computed  at  point  PA  for 
object  A  using  the  level  set  field.  The  normal  and  tangential  velocity  components  at  point 
PA  are  related  to  the  x-  and  y-direction  velocities  as: 
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v„=upnx+vpny  (46a) 

v,  =  upny  +  vpnx  (46b) 

Since  normal  velocities  are  continuous  across  the  interface  the  value  of  the  normal 
velocity  component  in  Eq.  (46a)  is  specified  to  be  that  of  the  object  B  at  that  point.  The 
tangential  velocity  component  in  Eq.  (46b)  is  assigned  to  be  the  extrapolated  value  from 
the  material  A.  With  these  values  of  vn  and  vt,  the  Cartesian  velocity  components  at  the 
point  Pa  are  then  obtained  from: 

up=v„nx+v,ny  (47a) 

vP=v„ny-vtnx  (47b) 

With  the  extrapolated  values  of  p  and  E,  pressure  is  obtained  from  equation  of  state. 

The  tangential,  normal,  and  shear  stresses  at  the  ghost  point  PA  are  computed  such  that 
the  continuity  of  normal  stress  is  enforced  at  the  interface.  The  components  of  stress 
oriented  with  the  local  tangent  and  normal  directions  are  computed  from: 

<ytt  =  «l,ASJV,A  +  ny,ASxx,A  ~  2nx,Any.AS xy,A  ~  P ,A  (48a) 

ann  =  nl,ASxx,B  +  nl,ASyy,B  +  lnx.Any.ASxy.B  -  P,B  (48b) 

=  nx<Any  A  (Syy  D  -  s„<B )  +  (nxA  -  nyA  )sxy  B  (48c) 

where  subscripts  t  and  n  denote  tangent  and  normal  directions  to  the  interface  at  point  O 
and  subscripts  A  and  B  indicate  the  material  from  which  the  data  are  obtained.  Note  that 
in  Eq.  (48a),  the  tangential  stress  is  computed  using  the  extrapolated  values  of  object  A  at 
PA.  This  implies  that  the  tangential  stress  remains  discontinuous  across  the  interface.  In 
Eqs  (48b)  and  (48c),  the  normal  and  shear  component  stresses  are  computed  using  field 
values  of  object  B  at  PA.  This  results  in  continuity  of  normal  and  shear  stresses  across 
the  interface. 

Converting  to  Cartesian  coordinates,  using  Eqs  (48a),  (48b),  and  (48c),  the  total  stresses 
at  the  ghost  point  PA  of  object  A  are  obtained: 


=  n2x,A(Tnn  +nl,A°tl-2nX'AnyACTnl 

(49a) 

=  n2yM<Tm+nlAatt^2nxMny:AaHt 

(49b) 

=  nx,Any,A^m  ~  (nx,A  ~  ny,A  )°»« 

(49c) 

Assuming  continuity  of  the  temperature  at  the  interface  between  the  materials,  the  value 
of  temperature  at  PA  is  assigned  the  value  of  temperature  in  object  B  at  the  grid  point 
coincident  with  PA.  The  value  of  equivalent  plastic  strain  is  extrapolated  from  the  interior 
of  object  A. 

When  interface  conditions  are  obtained  as  above  at  impact  boundaries  “overheating”  of 
the  material  can  result,  as  pointed  out  by  Glaister  (1988).  The  problem  arises  due  to  the 
collision  of  the  impinging  shock  on  a  solid  wall  and  the  reflected  shock  coming  off  a 
solid  wall  (Menikoff,  1994).  The  overshoots  in  the  density  and  temperature  are  numerical 
artifacts  that  arise  due  to  the  unphysical  dissipation  inherent  in  the  numerical  schemes 
and  the  inability  of  the  Euler  equations  to  conduct  heat  away  from  the  overheated  region 
into  the  wall  thus  causing  a  buildup  of  temperature  (Donat  and  Marquina,1996).  Pressure 
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and  velocity  appear  to  equilibrate  quickly  but  internal  energy  (temperature)  does  not  and 
the  use  of  the  e.o.s.  renders  the  density  value  inaccurate  also.  This  problem  has  been 
shown  in  a  previous  paper  (Udaykumar  et  al.,  2003)  to  exist  in  both  Eulerian  and 
Lagrangian  computations  of  impact.  A  fix  for  this  problem  has  been  suggested  by 
Fedkiw  et  al.  (1999b)  for  the  case  of  gases  and  other  materials.  In  the  present 
computations  with  solids  this  fix  has  not  been  applied  since  further  examination  of  this 
condition  is  required  to  determine  its  suitability  to  materials  with  strength,  since  the 
pressure  in  the  present  case  is  also  related  to  the  stresses  in  the  material. 


Type  2:  Material- void  interface 

For  a  material-void  interface,  the  physically  imposed  conditions  on  the  interface  are  that 
the  surface  tractions  be  negligible.  Therefore: 

n  •  (-pi  +  T)  =  0  (50) 

where  I  is  the  unit  tensor  and  T  is  the  deviatoric  stress  tensor.  In  the  ID  case  the  zero 
traction  condition  reduces  to: 


a*  =(\~P)±  =0  (51) 

This  condition  is  easily  applied  at  the  material-void  interface  in  ID  for  the  independent 
variable  sx  at  the  interface,  since  the  pressure  is  given  by  the  equation  of  state. 

In  the  two-dimensional  case,  implementation  of  this  boundary  condition  requires  care. 
To  apply  the  zero-traction  condition,  we  first  rotate  the  stress  tensor  as  follows.  Let 

a  =  AtgA  (52) 


where  6  = 


and  A  = 


nx  -ny 


are  the  rotated  stress  and  orientation  matrices 


J  L  ny  n*  J 

due  to  the  transformation  from  x-y  to  t-n  coordinates,  the  latter  coordinates  having  axes 
oriented  along  the  tangent  and  normal  to  the  interface.  The  physical  conditions  at  the 
interface  require  that  the  normal  and  shear  components  of  traction  be  zero,  while  the 
tangential  component  may  remain  discontinuous.  As  illustrated  in  Figure  4(b),  we  first 
compute  the  unit  normal  vector  at  the  ghost  point  P  (which  lies  in  the  void)  from  the 
level  set  field.  The  reflection  of  point  P  across  the  interface  lies  inside  the  material  at 
point  I.  Expanding  Eq.  (52),  the  stresses  oriented  along  the  normal  and  tangent  at  P  are 
obtained  at  the  point  I  by  interpolation  from: 


S’™,/  =  »ls„j  +"lsyyJ  +  2nxnysxy,i  ~  Pj 


(53a) 


~  (nx  ~  ny  )Sxy,l  +  nxHy  ( Syy,l  S xx.I )  (53b) 

Since  the  values  at  the  interface  for  all  tractions  are  zero  for  om  and  cr„, ,  linear 
extrapolation  is  performed  to  obtain  the  onn  and  anl  at  point  P,  using  the  values  at  I 
given  by  Eqs.  (53)  and  the  interfacial  values.  If  the  distance  between  point  E  and  the 
interface  is  less  than  a  certain  tolerance  (here,  we  use  O.lx  min(Ax,  Ay) ),  then  cr„„  and 

ont  are  then  taken  to  be  zero  at  the  point  P.  If  not,  then: 
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(54a) 


I 

® nn,P  ~  ^  ^ nnj 

l 

&nt,P  ~  d  Gnt,l 


(54b) 


The  stress  component  ou  is  not  required  to  be  continuous  across  the  interface  and  is 
obtained  using  the  extrapolated  ghost  values  at  P: 

&tt,P  ~  nxSyy,P  nysxx,p  ~  2 nxnySxy,P  ~  P,P  (^) 

Converting  to  Cartesian  components,  from  Eqs  (54a),  (54b),  and  (55),  the  stresses  at  the 
ghost  point  P  are  obtained  from: 

°xx,p  =  nl°nn,p  ~ 2 nxnyant  P  +  n2yattP  (56a) 

ayy,P  =  nlann,P  +  2nxny°n<,P  +  *1° 'tt.P  (56b) 

<?xy,P  =  nxny{Gm  P  -  o ttP )  +  (n2x  -  n2y)anl  P  (56c) 

All  other  variables  are  extrapolated.  Note  that  no  boundary  conditions  are  required  on 
the  void  side  of  the  interface  since  the  equations  are  not  solved  in  the  void  region. 


4.  Interface  Capturing  Scheme 

Discontinuities  in  properties  will  typically  exist  across  interfaces  between  interacting 
materials  and  will  need  to  be  maintained  during  their  interactions.  For  example,  the 
density  and  material  strength  properties  across  an  interface  could  be  different  by  several 
orders  of  magnitude  (e.g.  at  an  air-copper  boundary).  Traditional  purely  Eulerian 
numerical  methods  solve  the  governing  equations  by  adopting  the  single  domain 
approach,  i.e.  the  same  governing  equations  are  solved  at  each  node  in  the  domain, 
whether  void  or  material  node.  The  embedded  interfaces  and  accompanying 
discontinuities  are  then  accounted  for  by  using  different  techniques,  including  volume¬ 
averaging,  mixture  theories  (Benson,  1997),  or  by  using  “numerical”  delta  or  heaviside 
functions  (Peskin,  1977,  Sussman  et  al.,  1998).  As  an  alternative  to  such  Eulerian 
(diffuse  interface)  methods  it  is  possible  to  formulate  a  sharp-interface  technique  on  fixed 
meshes.  In  sharp-interface  approaches,  accurately  maintaining  the  physically  expected 
sharp  interfaces  between  materials  requires  sophisticated  interface  tracking  methods  and 
discretization  schemes  that  account  for  the  embedded  sharp  interfaces.  Such  approaches 
have  been  advanced  recently  by  various  researchers  (Fedkiw  et  al.,  1999a,  Leveque  et  al., 
1994,  Udaykumar  et  al.,  2003,  Tran  et  al.,  1999). 

In  previous  work  (Udaykumar  et  al.,  2003),  the  interface  was  tracked  explicitly  using 
marker  particles  and  parametrized  curves.  Such  an  approach  can  become  complicated  in  a 
3-dimensional  setting  or  when  the  interfaces  undergo  topological  changes.  An  alternative 
is  to  track  the  interfaces  implicitly  over  the  mesh  using  the  level  set  approach.  The 
advantage  of  the  level  set  approach  is  that  while  a  purely  Eulerian  update  is  employed  to 
advance  the  interface  via  the  level  set  (distance  function)  representation,  the  exact 
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location  of  the  interface  can  be  deduced  from  the  level  set  field.  Thus  the  level  set 
method  is  compatible  with  a  sharp-interface  treatment.  Furthermore,  the  level  set 
approach  easily  extends  to  3D  and  can  naturally  handle  merger/fragmentation  of  the 
interface.  However,  grid-supported  level  set  methods  build  in  an  entropy  condition  that 
leads  to  smoothing  of  sharp  comers  (cusps)  on  the  interface.  While  in  most  physical 
problems  this  feature  of  level  set  advection  can  be  advantageous,  as  unphysical  re-entrant 
interfaces  are  naturally  avoided,  in  the  cases  of  interest  to  the  present  work  this  aspect  of 
the  level  set  method  proves  to  be  a  hindrance.  This  is  because  we  wish  to  follow  the 
dynamics  of  solid  objects,  which  may  possess  comers  and  other  geometric  features  which 
must  be  carried  without  deterioration  over  the  mesh.  Recently,  a  Hybrid  Particle  Level 
Set  (HPLS)  method  was  presented  by  Enright  et  al.  (2001),  which  allows  transport  of 
comers  and  other  sharp  geometric  features  without  dissipation.  This  method  is  briefly 
described  in  the  following  section  and  the  results  from  the  calculations  using  the  method 
demonstrate  the  usefulness  and  appropriateness  of  the  method  for  the  class  of  problems 
solved  in  the  present  study. 


4.1  Level  Set  Method 


In  the  present  work  a  local-level  set  formulation  is  used  to  capture  the  moving  boundaries 
on  a  fixed  Cartesian  computational  grid.  The  method  introduces  a  continuous  scalar 
function  ®(3t,r) ,  where  the  position  of  the  actual  interface  is  identified  by  an  iso-value  of 
the  field,  i.e.  0  =  0.  If  multiple  objects  are  present,  a  separate  level  set  function  is  used 
to  describe  each  object. 


The  motion  of  the  interface  is  determined  by  a  velocity  field  u ,  which  can  depend  on 
position,  time,  and  geometry  of  the  interface.  The  motion  of  the  interface  can  thus  be 
described  by  a  scalar  convection  equation: 

—  +  «VO  =  0  (57) 

dt 

Since  we  are  only  interested  in  the  location  of  the  interface  ®(x,  t)  =  0,  the  above 
equation  only  needs  to  be  solved  locally  near  the  interface  (Adalsteinson  and  Sethian, 
1995;  Peng  et  al.,  1999).  Eq.  (57)  is  integrated  using  fourth-order  ENO  scheme  in  space 
and  a  third-order  Runge-Kutta  scheme  in  time. 

The  unit  normal  on  the  interface,  drawn  from  inside  the  material  (®  <  0)  to  outside  of 
the  material  (®  >  0) ,  and  the  curvature  of  the  interface  can  easily  be  expressed  in  terms 


of  ®(x,0: 


The  standard  reinitialization 
to  steady  state  the  equation: 


_  V® 
n  = 


and  k  =  V  • 

( v®"| 

0=0 

iNJ 

(58) 


0=0 


algorithm  maintains  the  signed  distance  property  by  solving 
<j>r  +  w  •  =  sign(0)  (59) 


where  r  is  the  fictitious  time, 
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sign(  0)  =  -=£=  and  w  =  5ign(<D)i^7  (60) 

fO  +  (Ax)2  lV^l 

with  the  initial  condition: 

#5,0)  =  0(5)  (61) 

However,  due  to  the  built-in  dissipation  in  the  ENO  scheme  and  the  reinitialization 
procedure,  excessive  regularization  occurs  at  under-resolved  regions.  Thus  sharp  comers 
are  rounded  off  and  objects  may  lose  mass  during  large  deformations.  Introducing 
Lagrangian  particles  in  combination  with  the  grid-based  level  set  can  remedy  this 
problem  as  shown  by  Enright  et  al.(2001). 

4.2  Hybrid  Particle  Level  Set  (HPLS)  Method 

The  hybrid  particle  level  set  method  (Enright  et  al.,  2001)  is  the  combination  of  an 
Eulerian  level  set  method  and  particle-based  Lagrangian  scheme.  The  particles,  of 
positive  and  negative  type,  are  placed  randomly  near  the  interface,  attracted  to  the  correct 
side  to  the  interface  (positive  for  O  >  0  and  negative  for  O  <  0)  and  passively  advected 
with  the  flow  using  the  equation: 

^  (62) 

at 

where  xp  is  the  position  of  the  particle,  which  is  the  center  of  a  sphere  of  radius  rp  and 
u(xp )  is  its  velocity.  Since  this  advection  procedure  entails  no  dissipation,  the  distance 

function  values  carried  by  the  Lagrangian  particles  can  be  considered  to  be  “correct”. 
Thus,  if  the  grid-based  level  set  function  differs  in  some  region  of  the  interface  from  the 
level  set  values  carried  by  the  particles,  the  grid-based  level  set  function  is  taken  to  be  in 
need  of  correction. 

While  the  particles  do  not  have  mass,  they  have  volume.  The  radius  of  each  particle  is 
bounded  by  minimum  and  maximum  values  based  upon  grid  spacing.  The  maximum  and 
minimum  radii  which  work  well,  according  to  Enright  et.  al.,  are: 

r^n  =0.1  min(Ax,  Ay,  A z)  (63) 

rn*x  =  °-5  min(Ax,  Ay,  A z)  (64) 

Initially  particles  of  both  signs  are  randomly  placed  within  a  band  of  3  x  max(Ax,  Ay,  Az) 
of  the  interface.  Particles  are  then  attracted  to  the  correct  side  of  the  interface  (i.e. 
positive  particles  to  the  O  >  0  side  and  negative  particles  to  the  <X>  <  0  side).  Finally,  each 
particle  radius  is  set  according  to: 

^max  ifSp<^Xp)>rimx 

rp  =  <  spO(xp )  if  rmin  <  sp<S>(xp )  <  r.  (65) 

/•min  if  SpO(xp)<rnin 

The  marker  particles  and  the  level  set  function  are  separately  integrated  forward  in  time. 
The  particles  that  escaped  to  the  wrong  side  of  the  interface  are  used  to  locate  possible 
errors  in  the  level  set  function  since  the  escaped  particles  indicate  that  the  characteristics 
have  most  likely  been  incorrectly  merged  through  regularization,  i.e.  the  level  set  has 
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computed  an  incorrect  weak  solution  (see  Figure  5).  Once  errors  are  located,  reduction  of 
error  is  carried  out  to  rebuild  the  grid-based  level  set  at  these  local  regions.  Detailed 
implementation  is  discussed  in  the  paper  by  Enright  et  al  (2001). 


The  error  correction  procedure  due  to  Enright  et  al.  is  as  follows.  Given  a  grid-based 
level  set  <X>  identify  the  set  of  escaped  positive  particles  E*  (a  escaped  positive  particle 
is  defined  as  a  positive  type  particle  that  happens  to  be  located  on  the  negative  side  of  the 
grid  level  set  as  illustrated  in  Figure  5).  Initialize  d>+  with  <t> ,  and  then  calculate, 

0+  =  max(0  ,C>+)  (66) 

VpeE* 


where  Op(3c)  is  defined  as: 

0,(x)  =  5,(r„-|*-.x,|)  (67) 

and  sp  is  the  sign  of  the  particle. 


Similarly,  to  calculate  a  reduced  error  presentation  of  the  O  <  0  region,  initialize  O 
with  ® ,  then  calculate, 

O'  =  max  (O  ,<!>')  (68) 

V/?e£" 


0>+  and  O'  are  then  merged  back  into  a  single  level  set  by  setting  O  equal  to  the  value 
of  0+  or  O"  which  is  least  in  magnitude  at  each  grid  point, 


i/|o+|  <|0' 

z/|o+|  >  lO' 


(69) 


The  minimum  magnitude  is  used  to  reconstruct  the  interface,  since  it  gives  priority  to 
values  that  are  closer  to  the  interface. 


The  above  procedure  for  error  correction  was  found  to  yield  somewhat  unsatisfactory 
results  for  certain  types  of  rigid  body  motions  when  the  meshes  used  were  not  sufficiently 
fine.  When  several  field  equations,  such  as  in  the  present  case,  need  to  be  solved 
simultaneously  with  the  level-set  equation,  the  grid  sizes  required  to  robustly  track  the 
sharp  comers  become  rather  prohibitive.  A  modification  to  the  error  correction  procedure 
of  Enright  et.  al.  was  therefore  made,  as  given  below,  which  appears  to  be  more  robust 
(although  somewhat  more  dissipative),  based  on  test  problems  run  using  modest  grid 
sizes. 

Given  a  level-set  function  O  and  a  set  of  escaped  positive  particles  E+ ,  we  initialize  the 
corrected  distance  function  values  at  the  4  surrounding  grid  points  0+  with  O ,  the  values 
obtained  by  advection  of  the  grid-based  level  set,  and  then  calculate: 

0+  =  max(|o.|,|o+|)  z/0+  >0  (70a) 

VpeE*  1  p'  1  1 

<D+  =-min(bl,b+|)  z/<D+  <0  (70b) 

Vpe£+  1  P''  1 
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For  all  the  escaped  particles  marked  as  positive  type  (Figure  5),  Eqs  70(a)  and  70(b) 
attempt  to  correct  the  interface  location  by  computing  new  level  set  values  at  the  four 
comers.  This  in  effect  moves  the  grid-based  interface  (i.e.  the  0=0  level)  so  that  the 
escaped  positive  particle  would  be  returned  to  the  correct  distance  on  the  positive  side  of 
the  interface.  Similarly,  a  correction  of  the  grid  values  can  be  applied  in  the  correction 
procedure  for  the  negative  escaped  particles.  Thus,  for  a  set  of  escaped  negative  particles 
E~ ,  we  initialize  O'  with  O  and  then  calculate 

O"  =  min(|0/,|,|o~|)  i/O'>0  (71a) 

O'  =  -max^O^O'])  i/O'<0  (71b) 

The  0+  and  O'  values  are  then  merged  back  to  obtain  the  distance  function  field  O  on 
the  grid  as  described  above. 

4.3  Detecting  and  Resolving  Collisions 

In  the  present  work,  material  interfaces  are  expected  to  collide  with  other  interfaces  or 
collapse  and  fragment.  Such  events  need  to  be  tracked  and  appropriate  interface 
conditions  applied  on  the  interacting  parts  of  the  interfaces.  In  the  current  framework, 
multiple  objects  with  different  material  properties  are  described  by  different  level  set 
functions.  Only  one  material  can  possess  a  given  computational  grid  point  at  any  time, 
whether  it  be  void  or  solid  material. 

The  collision  between  various  objects  is  determined  by  first  defining  the  properties  of 
each  level  set.  Each  level  set  (object)  is  associated  with  a  set  of  grid  points  which  straddle 
the  object  boundary,  i.e.  the  zero-contour  of  the  level-set.  In  2D,  these  grid  points  are 
identified  as  those  across  which  the  level  set  values  change  sign  from  negative  to 
positive.  Among  this  set,  the  grid  points  that  lie  within  the  material  defined  by  that 
particular  level  set  function  are  called  “boundary”  points  and  those  lying  just  outside  are 
classified  as  “ghost”  points.  This  set  of  points  is  stored  in  a  one-dimensional  array,  with 
mapping  assigned  to  access  the  grid  indexes.  At  any  given  time  step,  one  can  search 
through  this  set  of  grid  points  that  straddle  the  zero-level  for  each  object  to  determine 
whether  these  points  are  free  boundary  nodes  or  collision  nodes  as  described  below. 

Various  situations  can  arise  when  two  different  objects  move  toward  each  other  as  shown 
in  Figures  2(c,  i-iii).  The  distance  between  two  interfaces  is  easily  computed  using  the 
level  set  functions  associated  with  each  interface.  At  any  grid  point,  if  we  are  solving  for 
two  level  set  functions,  the  values  of  each  level  set  at  that  point  represent  the  distance  to 
their  respective  interfaces.  Therefore,  the  difference  in  these  values  represents  the 
distance  between  the  two  interfaces.  If  the  distance  between  two  approaching  level  sets 
is  less  than  a  specified  tolerance,  set  to  O.lx  min(Ax,Ay) ,  as  in  Figures  2(c)(ii  &  iii),  then 
the  ghost  point  is  marked  as  a  collision  boundary  point.  If  not  (Figure  2(c)(i)),  then  the 
ghost  point  is  marked  as  a  free  boundary  point.  This  process  is  repeated  for  each  level  set 
after  the  calculation  of  the  interface  velocity  for  each  level  set.  For  the  collision 
boundary  nodes,  each  level  set  interface  velocity  must  be  corrected  so  that  one  object  will 
not  penetrate  the  other.  For  two  level  sets,  each  defining  a  moving  boundary  under 
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impending  collision,  the  normal  velocity  at  the  interface  is  continuous  through  the  impact 
surface,  whereas  the  tangential  component  may  be  discontinuous  as  sliding  is  permitted. 
No  surface  friction  is  accounted  for  in  this  work  but  can  easily  be  included  in  the  present 
framework.  If  one  level  set  is  specified  as  moving  boundary  type  and  the  other  as  rigid 
stationary  boundary  type,  then  the  normal  component  of  velocity  is  zero  when  two  such 
interfaces  collide. 

Having  marked  the  ghost  points  as  collision  or  free  boundary  points,  the  boundary 
conditions  can  be  applied  at  the  interface  locations  as  described  in  section  3.3.  The  level 
set  velocities  are  then  obtained  and  the  interfaces  are  advected  to  their  new  positions. 
After  all  the  level  set  functions  are  evolved,  we  then  again  find  a  new  set  of  ghost  points 
and  determine  the  collision  nodes. 


5.  Results 

The  numerical  scheme  applied  to  evolve  the  flowfields  through  Eq.(14)  has  been  tested 
for  ID  problems  involving  impact  and  deformation  in  previous  work  (Udaykumar  et  al., 
2003).  Here  we  apply  the  scheme,  in  combination  with  the  particle-level-set  interface 
capturing  method  and  the  boundary  treatments  described  in  previous  sections  to  two- 
dimensional  problems.  The  results  are  compared  with  benchmark  solutions  and 
experimental  data. 

5.1  2D  Axisymmetric  Taylor  Impact  of  a  Copper  Rod 

In  the  2D  setting,  we  first  validate  the  methods  presented  above  by  testing  against  the 
well  known  Taylor  bar  benchmark  test  (Zhu  and  Cescotto,  1995).  A  cylindrical  rod  made 
of  copper  with  an  initial  radius  of  3.2  mm  and  a  length  of  32.4  mm  is  given  a  velocity  of 
227  m/s  (schematic  shown  in  Figure  6)  and  impacts  against  a  rigid  planar  surface.  The 
deformation  of  the  rod  is  presumed  to  be  axisymmetric.  The  rod  has  an  initial  density  of 
8930  kg/m3,  Young's  modulus  E  =  1  MGPa ,  Poisson’s  ratio  v  =  0.35 ,  and  yield  stress 
ay  =  400 MPa .  The  material  is  assumed  to  harden  linearly  with  a  plastic  modulus  of 

100 MPa .  The  calculations  are  carried  out  up  to  a  time  of  80/zs ,  at  which  point  nearly  all 
the  initial  kinetic  energy  has  been  dissipated  as  plastic  work  (Hallquist,  1987). 

Two  level  sets  are  initialized,  one  for  the  copper  rod  and  the  other  for  the  rigid  body. 
Since  the  bottom  rigid  body  is  not  moving,  there  is  no  need  to  initialize  it  with  particles. 
Therefore  particles  are  seeded  only  for  the  copper  rod.  Collision  between  two  bodies  is 
detected  using  the  algorithm  discussed  in  section  4.3,  and  the  calculations  are  carried  out 
with  the  material-material  and  material-void  boundary  conditions  applying  on  the 
collision  and  free  surfaces  of  the  rod,  as  discussed  in  section  3.4. 

A  grid  independence  study  is  performed  using  20x100,  40x200,  60x300,  and  80x400 
meshes.  Figure  7  shows  the  final  shapes  (at  80/zs)  for  all  the  meshes  mentioned.  It  is 
clear  that  the  solution  converges  with  grid  refinement  and  grid  independent  results  are 
approached  for  the  finest  mesh,  80x400. 
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Figure  8  shows  the  (grid-independent)  solutions  for  the  80x400  mesh  at  four  different 
times  after  impact,  i.e.  20  jjs ,  40  jus ,  60  jus ,  and  SO  jus .  It  can  be  observed  that  the  rod 
initially  extends  out  rapidly  at  the  point  of  impact  and  continues  to  extend  at  the  foot  up 
to  around  40 jus.  As  the  material  in  the  foot  hardens,  a  second  bulge  is  formed, 
extending  out  to  approximately  0.004  m  from  the  axis.  The  extent  of  the  bulge,  both  in 
the  radial  and  longitudinal  directions,  is  in  excellent  agreement  with  the  results  of 
Camacho  and  Ortiz  (1997).  The  rod  reaches  a  rest  state  around  the  time  of  SO  jus  in 
agreement  with  Camacho  and  Ortiz,  showing  that  the  conversion  of  kinetic  energy  to 
plastic  work  is  correctly  predicted. 

The  numerical  values  for  final  rod  length,  radius,  and  maximum  equivalent  plastic  strain 
are  compared  with  other  results  in  literature  in  Table  1.  Our  results  fall  squarely  in  the 
range  of  values  reported  in  literature.  The  slight  discrepancies  with  the  results  of 
Camacho  and  Ortiz  may  be  due  to  differences  in  the  treatment  of  the  rate-dependent 
plasticity  model  and  due  to  the  comparatively  coarser  mesh  in  the  present  case  near  the 
impact  surface  due  to  the  use  of  a  uniform  Cartesian  mesh.  At  the  moment  of  impact, 
extreme  stress  and  velocity  gradients  are  developed  in  the  rod  and  capturing  these 
features  may  require  a  very  fine  mesh.  Camacho  et  al.  used  a  moving  finite  element 
mesh  with  mesh  clustering  in  the  high-gradient  regions.  In  future  extensions  of  the 
present  work,  adaptive  local  refinement  techniques  will  be  implemented  in  the  current 
Cartesian  grid  framework  to  better  resolve  such  transient  features. 

The  evolution  of  velocity  magnitudes  along  with  velocity  vectors  are  shown  in  Figure 
9(a-d)  for  four  different  instants  after  impact.  Initially,  the  base  of  the  rod  expands  out 
rapidly,  forming  the  foot-like  shape.  After  around  40 jus ,  the  base  hardens  and  remains 
practically  stationary,  and  the  remaining  kinetic  energy  is  dissipated  in  the  work  required 
to  form  the  second  bulge,  above  the  foot  of  the  rod.  At  SO  jus ,  the  copper  rod  is  nearly  at 
rest. 

Figures  10(a-d)  shows  the  equivalent  plastic  strain  contours  at  four  different  times,  20 jus , 
40/zs ,  60 jjs  ,  and  80 jjs  .  One  can  see  the  development  of  large  plastic  strains  in  the  foot 
region,  and  a  local  minimum  in  the  plastic  strain  field  at  the  symmetry  axis,  centered 
between  3  and  4  mm  above  the  rigid  surface  at  the  symmetry  axis.  A  large  region  of  the 
rod  registers  negligible  plastic  strain.  These  features  are  clearly  shown  in  Figure  10(d). 
The  values  and  distribution  of  the  plastic  strain  are  in  very  good  agreement  with  the 
results  of  Camacho  and  Ortiz  (1997). 

The  maximum  temperature  of  around  500K,  is  achieved  close  to  the  base  at  around  the 
40 jjs  instant,  i.e.  when  the  plastic  deformation  of  the  base  reaches  a  maximum,  as 
expected.  This  is  the  point  in  time  where  nearly  all  the  plastic  work  done  in  deforming 
the  foot  has  been  converted  to  heat  and  raises  the  temperature  of  the  foot.  As  the  base  of 
the  foot  hardens,  the  maximum  rate  of  plastic  work  conversion  moves  up  to  form  the 
second  bulge  and  the  heat  is  dissipated  at  that  location.  However,  as  seen  from  Figure  10 
the  plastic  strain  accumulated  at  the  second  bulge  is  relatively  small. 
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5.2  2D  Axisymmetric  Penetration  of  Steel  Target  by  WHA  Long  Rod 

The  validation  of  our  method  for  two  deformable  objects  with  different  material 
properties  (a  case  requiring  2  different  level  sets)  is  carried  out  using  a  slender  tungsten 
heavy  alloy  (WHA)  rod  projectile  penetrating  an  initially  planar  target  made  of  a  steel 
plate  with  a  velocity  of  1250  m/s.  A  schematic  defining  the  problem  is  shown  in  Figure 
11.  Table  2  shows  the  material  properties  used  for  WHA  and  steel.  A  Johnson-Cook 
material  model  is  used  and  the  corresponding  strength  parameters  for  both  materials  are 
shown  in  Table  3.  Note  that  friction  between  the  two  impacting  surfaces  is  neglected  in 
these  calculations. 

Due  to  the  use  of  the  uniform  Cartesian  grid,  in  order  to  place  sufficient  mesh  points  to 
resolve  the  tungsten  rod,  a  very  large  mesh  size  is  required.  This  is  due  to  the  large  aspect 
ratio  of  the  rod  {length  I  radius  -  25).  Thus,  improving  the  resolution  in  the  radial 
direction  necessitates  the  use  of  fine  meshes  in  the  vertical  direction  also.  This  limitation 
can  be  overcome  by  use  of  adaptive  meshes  or  rectangular  Cartesian  meshes.  However,  if 
the  latter  is  used,  issues  concerning  inaccuracies  in  capturing  large  gradients  oriented 
arbitrarily  with  respect  to  a  stretched  mesh  will  need  to  be  addressed.  We  leave  these 
aspects  to  future  work.  To  enable  reasonable  computational  effort,  we  limit  the  axial 
dimension  of  the  problem  to  r  =  0.0125  and  study  what  effects  this  smaller  domain  has  on 
the  solutions. 

A  grid  independence  study  is  carried  out  using  3  meshes:  50x344,  75x516,  and  100x688. 
The  final  shapes  (at  80 /js  )  of  the  two  deformed  bodies  are  shown  in  Figures  12(a-c).  It  is 
seen  that  as  the  grid  is  refined  the  penetration  depth  as  well  as  the  length  of  the  upwelling 
WHA  material  (hereafter  called  ejecta)  become  larger.  At  the  finest  mesh,  100x688,  the 
solution  was  found  to  be  nearly  grid  independent.  Based  on  the  observed  trends,  further 
refinement  of  the  mesh  will  result  in  somewhat  larger  extents  of  the  ejecta  of  the  WHA 
material.  For  now,  we  will  use  the  100x688  mesh  results,  which  suffice  to  compare  the 
penetration  rate  and  depth  results  with  available  experimental  and  numerical  results 
(Anderson  et  al.,  1995  and  Camacho  and  Ortiz,  1997). 

Figures  12(d)  show  the  final  shape  of  the  two  deformed  bodies  at  80/zs  for  the  100x688 
mesh.  In  Figure  12(d),  we  see  that  the  right  side  of  the  steel  plate,  top  and  bottom,  barely 
moves  from  the  original  position.  In  reality  the  plate  used  in  experiment  has  a  larger 
radius  than  that  shown  in  Figures  12.  We  therefore  examined  whether  the  solutions 
obtained  from  our  calculations  on  the  truncated  domain  is  valid  by  extending  the  width  of 
steel  plate  from  0.0125  m  to  0.02  m,  while  maintaining  the  same  grid  density  as  the 
100x688  mesh.  This  extended  domain  case  for  the  steel  plate  (r  =  0.02  m)  is  shown  in 
Figures  13(a-d)  with  evolution  of  the  two  deformed  bodies  at  times  of  20 /js  ,  40 jus , 
60 /is ,  and  80/zs  for  a  mesh  of  160x688  points.  Note  that,  due  to  the  large  aspect  ratio  of 
the  problem,  even  for  this  fairly  large  mesh  size,  there  are  only  16  mesh  points  within  the 
initial  rod  radius.  Nevertheless,  the  results  obtained  agree  very  well  with  the  benchmarks 
and  with  experiments  as  demonstrated  below.  The  top  and  bottom  surfaces  of  the  steel 
plate  do  not  deform  much  at  r  >  0.01  m.  This  behavior  agrees  well  with  calculations 
reported  by  Camacho  and  Ortiz.  The  ejected  length  of  WHA  material  (Figure  13(d))  is 
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slightly  shorter  compared  to  that  for  narrower  domain  (r  =  0.0125  m)  in  Figure  12(d). 
This  is  perhaps  due  to  the  absorption  and  dissipation  of  energy  within  the  larger  mass  of 
material  in  the  plate  in  the  extended  domain  case.  In  particular,  the  top  surface  of  the  steel 
plate  appears  to  deform  more  for  the  extended  domain  case  shown  in  Figure  12(d). 
However,  the  main  features  of  the  computed  solutions  appear  to  be  relatively  insensitive 
to  the  domain  size  and  thus  the  truncation  of  the  domain  to  facilitate  computational 
economy  does  not  appear  to  be  very  significant. 

The  evolution  of  equivalent  plastic  strain  is  shown  in  Figures  14(a-d)  for  the  extended 
domain  with  a  160x688  mesh.  The  maximum  equivalent  plastic  strain  is  found  to  be 
around  4.5,  occurring  mostly  near  the  impact  surfaces.  The  values  of  equivalent  plastic 
strain  are  higher  in  the  WHA  material  compared  to  those  in  the  steel  material.  The 
plastic  strains  obtained  by  Camacho  and  Ortiz  using  Lagrangian  finite  element  method 
with  an  adaptive  mesh  agree  very  well  with  the  present  results,  both  in  terms  of  the 
magnitude  and  distributions  of  the  plastic  strains.  In  particular,  a  trough  in  the  plastic 
strain  distribution  is  noticed  in  both  our  results  as  well  as  those  of  Camacho  and  Ortiz  and 
occurs  near  the  bottom  surface  in  the  steel  plate  at  the  symmetry  axis,  as  seen  in  Figure 
14(e).  The  ejection  length  of  the  WHA  material  is  higher  in  the  Camacho  and  Ortiz 
calculations  when  compared  to  our  results.  However,  the  resolution  of  the  ejected  region 
afforded  by  the  mesh  used  in  the  present  calculations  is  too  low,  with  just  3  mesh  points 
across  the  vertically  oriented  trails  of  the  ejecta.  The  grid  refinement  study  performed 
above  indicates  that  as  the  mesh  is  refined  further  the  length  of  the  ejecta  will  increase. 
As  shown  below,  at  the  current  mesh  resolution,  the  overall  penetration  characteristics 
and  material  deformation  are  adequately  predicted. 

Figures  15(a-d)  shows  the  v-component  velocity  contour  at  different  times,  20 /js ,  40 /is  , 
(,0/js ,  and  %0/js,  for  the  extended  domain  with  a  160x688  mesh.  The  maximum  positive 
v-component  velocity  is  observed  around  40 /js  ,  occurring  in  the  ejecting  mass  of  the 
WHA  material.  Around  80 /js,  the  rod  comes  to  rest  and  only  small  residual  velocities 
remain. 

The  maximum  temperatures  occur  around  the  impacted  surfaces  as  this  the  region  of 
maximum  rate  of  conversion  of  plastic  work  to  heat.  The  recorded  maximum 
temperature  is  around  1300  K  in  the  WHA  material,  below  the  melting  temperature  of 
1777  K  for  WHA  and  1723  K  for  steel.  The  largest  temperature  occurs  at  around  40 /js  , 
and  decreases  as  the  rod  goes  to  rest  state.  This  shows  that  the  largest  plastic  work  done 
occurs  before  this  time  of  40 jus . 

Figure  16  shows  the  projectile  nose  and  tail  trajectories  as  a  function  of  time,  for  the 
extended  domain  case,  and  is  compared  with  the  superposed  results  from  experiment 
(Anderson  et  al.,  1995)  and  from  Camacho  and  Ortiz.  Also  plotted  are  the  original  rear 
and  impact  surfaces.  Our  results  show  reasonable  agreement  with  those  of  experiment 
and  Camacho  and  Ortiz.  The  tail  trajectory  is  in  much  better  agreement  as  its  surface 
experiences  less  extreme  conditions  during  impact  and  penetration.  The  present 
calculation  predicts  the  penetration  depth  in  good  agreement  with  experiments.  Despite 
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the  marginal  resolution  of  the  ejected  trails,  the  overall  penetration  and  deformation 
behavior  is  predicted  in  good  accord  with  the  adaptive  finite  element  simulations  of 
Camacho  and  Ortiz  (1997). 

Figure  17  shows  the  projectile  nose  (upper  curves)  and  tail  (lower  curves)  velocity 
histories,  for  the  extended  domain  case,  compared  against  the  results  from  experiment 
(Anderson  et  al.,  1995)  and  from  Camacho  and  Ortiz.  Note  that  since  the  results  from 
Camacho  and  Ortiz  display  considerable  noise  in  the  initial  phases  after  impact,  the 
plotted  values  were  hand-picked  to  show  only  a  few  representative  points  after  smoothing 
their  curves.  Similar  conclusions  can  be  drawn  as  in  Figure  16,  the  velocity  data  obtained 
from  the  present  calculations  being  in  good  agreement  with  the  benchmarks  with  the  tail 
results  being  in  better  agreement  than  nose  results.  The  oscillatory  nature  of  the  results, 
seen  in  both  the  present  and  the  FEM  calculations  of  Camacho  and  Ortiz  is  perhaps  due 
to  the  transient  wave  phenomena  in  the  initial  stages  of  impact. 

5.3  2D  Axisymmetric  Void  Collapse  in  a  Copper  Matrix 

In  this  test  case,  a  spherical  void  with  a  radius  of  1  micrometer  within  a  copper  matrix 
undergoes  axisymmetric  deformation  due  to  a  propagating  shock  created  by  imposing  a 
particle  velocity  at  the  bottom  boundary  (i.e.  simulating  a  piston  moving  at  specified 
speeds).  A  schematic  of  the  problem  is  shown  in  Figure  18.  The  impulsive  velocity 
supplied  at  the  bottom  boundary  seeks  to  simulate  the  effect  of  impact  at  the  bottom 
surface  of  the  copper  plate.  The  material  response  is  modeled  using  a  Johnson-Cook 
model,  with  constants  applicable  to  copper  (Table  3).  The  wave  moves  through  the 
copper  matrix  and  is  transmitted  out  through  the  upper  boundary.  The  right  boundary  is 
subjected  to  symmetric  conditions,  whereas  the  left  boundary  is  the  axis  of  symmetry. 
The  simulation  is  carried  out  for  three  different  initial  piston  velocities:  50  m/s,  200  m/s, 
and  500  m/s,  using  a  200x400  mesh.  The  void  collapse  phenomenon  has  implications  for 
the  initiation  process  in  energetic  materials  (Bowden  and  Yoffe,  1952,  Carroll  and  Holt, 
1972,  Kang  et  al.,  1992,  Menikoff  and  Kober,  1999).  Here  we  examine  the  behavior  of 
the  void  in  copper,  for  which  the  material  parameters  are  well  characterized.  Void 
collapse  can  occur  in  different  modes  depending  on  the  strength  of  the  impinging  shock. 
As  the  shock  strength  increases  the  void  collapse  goes  from  a  nearly  spherical  (visco¬ 
plastic)  mode  to  a  jet  (hydrodynamic)  mode  where  the  lower  surface  of  the  void  forms  a 
jet  which  impacts  on  the  upper  surface  at  high  velocity.  The  criterion  for  the  transition 
from  the  visco-plastic  to  hydrodynamic  mode  is  provided  by  the  analysis  of  Khasainov 
(1981)  in  terms  of  the  ratio  of  shock  passage  time  to  the  void  deformation  time  scales. 
When  the  shock  passage  time  scale  is  larger  than  the  void  collapse  time  scale  the  mode  of 
collapse  is  visco-plastic;  when  the  shock  passage  time  scale  is  comparable  with  the  void 
collapse  time  scale  the  mode  of  collapse  is  hydrodynamic.  This  latter  mode  is 
characterized  by  the  formation  of  a  jet  of  material  which  issues  from  the  lower  side  of  the 
void  and  impacts  the  upper  side,  leading  to  large  temperatures  on  the  impact  location  due 
to  dissipation  of  kinetic  energy  of  the  jet.  This  transition  is  displayed  in  the  following 
results. 
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For  the  low  impact  velocity  case,  i.e.  50m/s,  Figures  19(a-d)  show  the  void  at  t  =  2.52e-9 
s,  after  which  fiirther  deformation  nearly  ceases.  The  figures  show  the  final  shape  (19a), 
pressure  contour  (19a),  temperature  contours  (19b),  and  velocity  magnitude  contour  (19b) 
along  with  streamlines.  The  maximum  temperature  achieved  is  around  404  K,  occurring 
around  the  void.  The  collapsing  process  is  not  concentric  as  the  original  void  center  was 
at  (0,2e-06)  m.  However,  the  void  does  not  markedly  deviate  from  its  original  spherical 
shape.  No  jetting  phenomenon  is  observed,  although  the  maximum  v-component  velocity 
(~  100  m/s)  is  at  the  lower  surface  at  the  axisymmetric  coordinate.  The  kinetic  energy 
developed  due  to  the  imposed  shock  is  dissipated  through  plastic  work.  Figure  20  shows 
the  evolution  of  shapes  of  the  void  at  equal  time  intervals.  Most  of  the  deformation 
occurs  on  the  lower  surface,  and  during  the  earlier  phase  of  the  process. 

Results  for  imposed  boundary  velocity  of  200  m/s,  i.e.  intermediate  shock  strength,  are 
shown  in  Figures  21(a-d)  at  t  =  9.45e-10  s.  This  void  collapses  completely  and  the  instant 
shown  is  shortly  before  complete  void  collapse.  The  figures  show  the  final  shape  (21a), 
pressure  contour  (21a),  temperature  contour  (21b),  and  velocity  magnitude  contour  (2 Id) 
along  with  streamlines.  It  is  observed  that  in  this  case,  jetting  occurs  to  some  extent,  with 
maximum  jetting  velocity  around  2600  m/s.  The  maximum  temperature  in  the  material 
during  the  collapse  is  around  1570  K,  which  is  above  the  melting  temperature  of  1400  K 
for  copper.  The  highest  temperatures  occur  around  the  collapsing  void  as  this  is  where 
highest  rate  of  plastic  work  is  localized. 

Figure  22  shows  the  evolution  of  the  shapes  during  the  void  collapse  process  plotted  at 
equal  time  intervals.  The  jetting  phenomenon  is  clearly  observed,  eventually  breaking  the 
void  at  the  symmetry  axis  into  two  smaller  voids,  which  eventually  disappear.  This 
simulation  demonstrates  the  ability  of  the  method  to  handle  extreme  topological  changes. 

For  an  imposed  particle  velocity  of  500  m/s,  i.e.  the  highest  shock  strength,  the  flow  field 
at  three  different  times  during  the  collapse  process,  i.e.  3.78  x  10',0s ,  5.67  x  10_l05 ,  and 
6.93  x  10'105 ,  are  shown  in  Figures  23, 24,  and  25.  Final  shapes  (a),  pressure  contour  (a), 
temperature  contour  (b),  and  velocity  magnitude  contours  (b)  along  with  streamlines  are 
plotted  in  each  figure.  For  this  high  strength  of  shock,  a  distinct  jet  can  be  observed,  with 
the  highest  velocity  of  the  jetting  material  being  around  4130  m/s,  much  higher  than  the 
imposed  particle  velocity,  and  the  maximum  temperature  recorded  during  the  collapse 
process  is  around  2110  K,  well  above  the  melting  temperature  of  copper.  One  can  see 
that  the  void  would  be  completely  collapsed  before  the  shock  wave  reaches  the  top 
boundary  of  the  domain,  i.e.  The  shock  passage  time  in  this  case  is  comparable  with  the 
void  collapse  time.  This  case  represents  the  hydrodynamic  mode  of  collapse  in 
agreement  with  the  criterion  mentioned  above.  The  maximum  temperatures  are 
concentrated  around  the  void,  and  increase  from  426  K  (at  3.78xlO“los)  to  2110  K  (at 
6.93  x  10_I° s ),  as  the  plastic  work  is  dissipated  as  heat. 

Figure  26  shows  the  evolution  of  the  void  collapse  process.  One  can  see  a  clear  jetting 
phenomenon  after  the  time  of  5.67  x  10'lo5 .  Following  the  final  shape  shown,  the  void  is 
completely  collapsed  and  thus  the  void  space  disappears.  The  computational  method  is 
capable  of  simulating  the  entire  process  of  void  deformation,  collapse  and  disappearance. 
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6.  Conclusions 

We  have  developed  a  numerical  technique  designed  to  solve  problems  involving  high¬ 
speed  multimaterial  interactions  occurring  in  impact,  penetration  and  collapse.  The 
distinguishing  feature  of  the  approach  presented  here  is  that  the  shortcomings  associated 
with  traditional  Eulerian  and  Lagrangian  methods  are  avoided.  The  methodology  can 
easily  be  extended  to  three-dimensions.  The  technique  can  handle  the  following  physical 
phenomena  typical  of  impact  problems: 

1.  Large  deformations,  including  fragmentation  and  merger  of  the  materials.  All 
boundaries  are  treated  in  a  sharp  manner.  The  grid  remains  fixed  however  and 
problems  associated  with  maintaining  a  high-quality  grid  under  large 
deformations  do  not  arise. 

2.  Nonlinear  wave-propagation  and  the  development  of  shocks  in  materials 
governed  by  rate-dependent  plasticity.  The  nonlinear  waves  are  tracked  using 
high-resolution  shock  capturing  schemes  implemented  on  the  fixed  grid. 

3.  Accurate  elasto-plastic  behavior  of  the  material  during  impact.  This  aspect  is 
captured  for  material  response  governed  by  the  Johnson-Cook  model  and  a  radial 
return  algorithm  to  maintain  consistency  with  the  yield  surface. 

For  2-dimensional  problems  the  hybrid  particle  level  set  method  is  used  to  track 
boundaries  with  sharp  comers  that  are  carried  without  deterioration  through  the  large 
deformations  of  the  materials.  Benchmark  calculations  for  the  multi-dimensional  case 
including  axisymmetric  Taylor  bar  impact  and  penetration  of  a  Tungsten  rod  into  steel 
plate  show  excellent  agreement  with  moving  finite  element  solutions.  Qualitative 
agreement  with  theory  is  shown  for  void  collapsing  process  in  an  impacted  material 
containing  a  spherical  void.  The  method  has  thus  been  shown  to  be  suitable  for 
applications  involving  high-velocity,  multimaterial  impacts  leading  to  large  strain-rates, 
nonlinear  elasto-plastic  waves  and  topological  changes. 
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Case 

Final  length  (mm) 

Max.  equivalent 
plastic  strain 

Current 

21.15 

7.15 

2.86 

Camacho  &  Ortiz 

21.42-21.44 

7.21-7.24 

2.97-3.25 

DYNA2D 

21.47 

7.13 

3.05 

Zhu  &  Cescotto 

21.26-21.40 

6.97-7.18 

2.75-3.03 

Table  1 .  Comparison  of  results  for  Taylor  bar  impact  problem. 
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Material 

P 

(kg/m3) 

E  (GPa) 

u 

C 

(W/m/K) 

k 

(J/kg/K) 

r 

Co  (m/s) 

Tm  (K) 

Tungsten 
heavy  alloy 

17600 

200 

0.29 

477 

38 

1.43 

4030 

1777 

High- 

strength  steel 

7850 

323 

0.30 

134 

75 

1.16 

4570 

1723 

Copper _ 

8930 

117 

0.35 

383.5 

401 

2.0 

3940 

1358 

Table  2.  Material  properties  and  Mie-Griineisen  equation  of  state  parameters. 


Material 

Y0  (GPa) 

B  (GPa) 

n 

C 

m 

G  (GPa) 

Tungsten 
heavy  allow 

1.51 

0.177 

0.12 

0.016 

1.00 

124.0 

High-hard 

steel 

1.50 

0.569 

0.22 

0.003 

1.17 

77.3 

Copper 

0.400 

0.100 

1.0 

0.025 

1.09 

43.33 

a,  =  (Y„  +  Bs‘r)[ 1  +  Cln(£,/4)](1  -  S'"),  4  =  1.0s'1, 


8-,T~T> 


T.-T0 


T„  =  294 K. 


Table  3.  Constitutive  parameters. 
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Immersed 

interface 


Figure  1  (a)  Illustration  of  grid  point  and  grid  face  definitions  for  discretization  of 
governing  equations.  H,+  and  FT  are  derivatives  of  the  interpolating  function  evaluated  from 
the  left  and  right  stencils  respectively,  (b)  Grid  point  and  grid  face  definitions  for  evaluation 
of  fluxes  in  the  presence  of  an  immersed  boundary,  (c)  Example  of  profile  of  q  variable 
along  x. 
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material-void  region 


material-material  region 


Figure  2.  a)  Illustration  of  boundary  types  when  multiple  bodies  interact,  b)  Zoom-in 
view  of  two  bodies  in  contact  show  ghost  and  collision  nodes,  c)  Collision  scenarios  of 
two  bodies  in  ID. 
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Neighbor  nodes  of  I2 


Figure  3.  Illustration  of  extrapolation  process. 


A 

nB 


Figure  4.  Illustration  of  boundary  conditions  for  a)  Material-Material 
boundaries  and  b)  Material- Void  boundary. 
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Escaped  negative  particles 


Figure  5.  Illustration  of  escaped  positive  and  negative  particles  used  in  error 
correction  procedure. 
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40 


3.2  mm 


Figure  6.  Schematic  of  the  Taylor  bar  impact  problem. 


41 


z  im> 
0.035 

0030 

0025 


0020 

0.015 

0.0L0 

0.005 


0000 

0.0000  00050 


J.  J  ]  I  (m) 
00100 


Figure  7.  Grid  refinement  study  in  Taylor  bar  impact  calculation,  with  a 

velocity  of  227  m/s.  Final  shapes  of  20x100  mesh  ( - ),  40x200 

mesh  ( - ),  60x300  mesh  ( - ),  and  80x400  mesh  ( - ). 
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Figure  8.  Shapes  of  Taylor  bar  impact  calculation  at  different 
times  for  80x400  mesh,  with  a  velocity  of  227  m/s.  At  20  /is 
( - ),  40  /is  ( - ),  60  /is  ( - ),  80  /is  ( - ). 
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Figure  10.  Velocity  magnitude  contour  along  with  velocity  vector  at  different  times 
using  80x400  mesh,  a)  20  [is,  b)  40  /us,  c)  60  jus,  and  d)  80  /us,  with  a  velocity  of  227 
m/s. 
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Diameter  =  0.4  cm 


-H  h-/ 


Figure  1 1 .  Schematic  of  a  tungsten  rod  penetrating  a  steel  plate  target. 


a)  b) 


c) 


Figure  12.  Grid  refinement  study  of  a  tungsten  rod  penetrating  a  steel  plate 
with  a  velocity  of  1250  m/s.  a)  50x344  mesh,  b)  75x516  mesh,  and  c) 
100x688  mesh. 
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Figure  13.  Sha 
with  a  velocity 


a)  b)  c) 


0000  0  010  0.020 


Figure  14.  Equivalent  plastic  strain  contour  of  a  tungsten  rod  penetrating  a  steel  plate 
calculation  using  160x688  mesh,  with  a  velocity  of  1250  m/s.  a)  20  [is,  b)  40  [is,  c)  60 
[is,  d)  80  [is. 
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a)  b) 


Figure  15.  V-component  velocity  contour  of  a  tungsten  rod  penetrating  a 
steel  plate  calculation  using  160x688  mesh,  with  a  velocity  of  1250  m/s. 
a)  20  iis,  b)  40  /ts,  c)  60  ns,  and  d)  80  ixs. 
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Figure  16.  Projectile  nose  and  tail  trajectories  of  a  tungsten  rod  penetrating  steel 
plate  using  a  160x688  mesh,  with  a  velocity  of  1250  m/s. 
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Figure  17.  Projectile  nose  and  tail  velocity  histories  of  a  tungsten  rod  penetrating 
steel  plate  using  a  160x688  mesh,  with  a  velocity  of  1250  m/s. 
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Figure  19.  Void  collapse  calculation  using  200x400  mesh  in  a  copper  matrix,  with  a  diameter 
=  1  fjm ,  V  =  50  m/s,  and  t  =  2.52x  10'9  s.  a)  left  -  shape,  right  -  pressure  contour,  b)  left  - 
temperature  contour,  right  -  velocity  magnitude  contour  and  velocity  streamline 
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Figure  20.  Shapes  at  different  times  of  void  collapse  process  using  200x400 
mesh  in  a  copper  matrix,  with  a  diameter  =  1  jjm  and  V  =  50  m/s. 
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Figure  21.  Void  collapse  calculation  using  200x400  mesh  in  a  copper  matrix,  with  a  diameter 
=  1  /m,V  =  200  m/s,  and  t  =  9.45  x  10~10  s.  a)  left  -  shape,  right  -  pressure  contour,  b)  left  - 
temperature  contour,  right  -  velocity  magnitude  contour  and  velocity  streamline. 
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Figure  22.  Shapes  at  different  times  of  void  collapse  process  using  200x400 
mesh  in  a  copper  matrix,  with  a  diameter  =  1  /mu  and  V  =  200  m/s. 
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Figure  23.  Void  collapse  calculation  using  200x400  mesh  in  a  copper  matrix,  with  a  diameter 
=  l  jum,W  =  500  m/s,  and  t  =  3.78x  10'10  s.  a)  left  -  shape,  right  -  pressure  contour,  b)  left  - 
temperature  contour,  right  -  velocity  magnitude  contour  and  velocity  streamline 
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Figure  24.  Void  collapse  calculation  using  200x400  mesh  in  a  copper  matrix,  with  a  diameter 
=  1  /jm ,  V  =  500  m/s,  and  t  =  5.67x  10"'°  s.  a)  left  -  shape,  right  -  pressure  contour,  b)  left  - 
temperature  contour,  right  -  velocity  magnitude  contour  and  velocity  streamline 
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Figure  25.  Void  collapse  calculation  using  200x400  mesh  in  a  copper  matrix,  with  a  diameter 
=  1  im ,  V  =  500  m/s,  and  t  =  6.93x10“'°  s.  a)  left  -  shape,  right  -  pressure  contour,  b)  left  - 
temperature  contour,  right  -  velocity  magnitude  contour  and  velocity  streamline 
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Figure  26.  Shapes  at  different  times  of  void  collapse  process  using  200x400 
mesh  in  a  copper  matrix,  with  a  diameter  =  1  fjm  and  V  =  500  m/s. 
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